#!/bin/bash
#从 kingdom 导出层位然后 line drawing,由于 horizon 需要新建很多层位,因此用断层表示层位,导出的都是断层格式(GeoQuest-CARD IMAGE 7),注意左对齐
function preset()
{
    #=================================
    fig_name=25_line_drawing
    fig_fmt=png
    #=================================
    gmt set MAP_FRAME_TYPE=plain
    lon_min=-72
    lon_max=-69.5
    lat_min=-22.2
    lat_max=-18.4
    R=$lon_min/$lon_max/$lat_min/$lat_max
    J=M5.5ih
}

preset

gmt begin $fig_name $fig_fmt E1000
# === (a) 25_line_drawing ===
origin_x_min=5300 #set min_cmp num 04
origin_x_max=25527 #set max_cmp num 04
new_x_min=`echo "scale=2; $origin_x_min * 6.25 / 1000 " | bc`
new_x_max=`echo "scale=2; $origin_x_max * 6.25 / 1000 " | bc`
origin_d_min=0
origin_d_max=40
xylimits_cmp=`echo $origin_x_min/$origin_x_max/$origin_d_min/$origin_d_max`
trench_cmp=7700
trench_km=`echo "scale=2; $trench_cmp * 6.25 / 1000 " | bc`
left_trench=`echo "scale=2; $new_x_min - $trench_km " | bc`
right_trench=`echo "scale=2; $new_x_max - $trench_km " | bc`
xylimits_km=`echo $left_trench/$right_trench/$origin_d_min/$origin_d_max` #6900*6.25/1000
height=1
vert_exa=1
width=`echo "scale=2; $height * (${origin_x_max} - ${origin_x_min}) * 6.25 *0.001  / ((${origin_d_max} - ${origin_d_min})* ${vert_exa})" | bc`
for v in 25  #读出起始终止坐标
do
    x1=$(awk '{if($1=='$origin_x_min') print $2}' ../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
    y1=$(awk '{if($1=='$origin_x_min') print $3}' ../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
    awk '{if($1=='$origin_x_min') print $2,$3}' ../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt > ${v}_start_utm.txt #utm
    x11=$(awk '{if($1=='$origin_x_min') print $2}' ../../gmt/P1_basemap_shade/MGL1610MC${v}.latlong) #latlong
    y11=$(awk '{if($1=='$origin_x_min') print $3}' ../../gmt/P1_basemap_shade/MGL1610MC${v}.latlong) #latlong
    x2=$(awk '{if($1=='$origin_x_max') print $2}' ../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
    y2=$(awk '{if($1=='$origin_x_max') print $3}' ../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
    awk '{if($1=='$origin_x_max') print $2,$3}' ../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt > ${v}_end_utm.txt #utm
    x21=$(awk '{if($1=='$origin_x_max') print $2}' ../../gmt/P1_basemap_shade/MGL1610MC${v}.latlong) #latlong
    y21=$(awk '{if($1=='$origin_x_max') print $3}' ../../gmt/P1_basemap_shade/MGL1610MC${v}.latlong) #latlong
done

for v in 25  #seafloor
    do
        gmt basemap -R$x1/$x2/$origin_d_min/$origin_d_max -JX${width}i/-${height}i -Bya10f2+l"Depth (km)" -BWeSn \
        --MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p
        awk '{print $1,$3/1000*-1}' ../slice_seismic/seafloor_val_${v}.utm > tmp.txt
        gmt plot tmp.txt -W0.1p,black  
    done
awk '$3~/[[:digit:]+]/{print $1,$4/1000,$3}' 25_basement.dat > tmp0.txt #basement
k=$(awk 'END{print $3}' tmp0.txt) #the max num of the layers
i=0
while [ ${i} -le ${k} ]; do
awk '$3 == '${i}'{print $1,$2}' tmp0.txt > tmp.txt
gmt plot tmp.txt -W0.1p,black
i=$(($i + 1))
done
awk '$3~/[[:digit:]+]/{print $1,$4/1000,$3}' 25_horizon.dat > tmp0.txt #horizon
k=$(awk 'END{print $3}' tmp0.txt) #the max num of the layers
i=0
while [ ${i} -le ${k} ]; do
awk '$3 == '${i}'{print $1,$2}' tmp0.txt > tmp.txt
gmt plot tmp.txt -W0.1p,magenta
i=$(($i + 1))
done
awk '$3~/[[:digit:]+]/{print $1,$4/1000,$3}' 25_f.dat > tmp0.txt #fault
k=$(awk 'END{print $3}' tmp0.txt) #the max num of the layers
i=0
while [ ${i} -le ${k} ]; do
awk '$3 == '${i}'{print $1,$2}' tmp0.txt > tmp.txt
gmt plot tmp.txt -W0.1p,red
i=$(($i + 1))
done
awk '$3~/[[:digit:]+]/{print $1,$4/1000,$3}' 25_rotat.dat > tmp0.txt #rotat
k=$(awk 'END{print $3}' tmp0.txt) #the max num of the layers
i=0
while [ ${i} -le ${k} ]; do
awk '$3 == '${i}'{print $1,$2}' tmp0.txt > tmp.txt
gmt plot tmp.txt -W0.1p,grey
i=$(($i + 1))
done
gmt psbasemap -R$xylimits_km -JX${width}i/-${height}i -Bxf5a10+l"Distance to trench (km)" -BS \
--MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p
gmt psbasemap -R$xylimits_cmp -JX${width}i/-${height}i -Bxf100a2000+l"CMP number" -BN \
--MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p
gmt text -F+cBL+f5p,1 -Dj0.1c/0.1c <<EOF
MC 25
EOF
gmt text -F+cBL+f5p,1 -Dj0.1c/0.3c <<EOF
VE=1
EOF
gmt text -F+cTL+f5p,1 -Dj0.1c/0.1c <<EOF
W
EOF
gmt text -F+cTR+f5p,1 -Dj0.1c/0.1c <<EOF
E
EOF

gmt end show
rm tmp*
rm gmt*